#include"phasefunction.h"
inline Float HenyeyGreenstein::phase(Float cosTheta)
{
	Float denom = 1 + g * g + 2 * g*cosTheta;
	return (1 - g * g)*Inv4Pi / (denom*sqrtf(denom));
}
HenyeyGreenstein::HenyeyGreenstein(Float g) : g(g) {}

Float HenyeyGreenstein::p(P3 &wo, P3 &wi)
{
	return phase(wo.dot(wi));
}

Float HenyeyGreenstein::sample_p(P3 &wo, P3 *wi, Float *rands)
{
	// Compute $\cos \theta$ for Henyey--Greenstein sample
	Float cosTheta;
	if (abs(g) < eps)
		cosTheta = 1 - 2 * rands[0];
	else {
		Float sqrTerm = (1 - g * g) / (1 - g + 2 * g * rands[0]);
		cosTheta = (1 + g * g - sqrTerm * sqrTerm) / (2 * g);
	}

	// Compute direction _wi_ for Henyey--Greenstein sample
	Float sinTheta = sqrtf(max((Float)0, 1 - cosTheta * cosTheta));
	Float phi = 2 * Pi * rands[1];
	Float cosPhi = cos(phi), sinPhi = sin(phi);
	P3 wi_local(sinPhi*cosTheta, sinPhi*sinTheta, cosPhi);

	//local coordinate system with wo being z-axis
	P3 u = ((fabs(wo.x) > 0.1f ? P3(0, 1, 0) : P3(1, 0, 0)).cross(wo));
	P3 v = wo.cross(u);
	*wi = u * wi_local.x + v * wi_local.y - wo * wi_local.z;
	return phase(-cosTheta);
}